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Abstract 

This paper presents a three-dimensional, energy based, anisotropic, stiffness reduction, progressive damage model 
for composite materials and composite material constituents. The model has been implemented as a user-defined 
constitutive model within the Abaqus finite element software package and applied to simulate the nonlinear behavior 
of a damaging epoxy matrix within a unidirectional composite material. Three different composite microstructures 
were considered as finite element repeating unit cells, with appropriate periodicity conditions applied at the bound- 
aries. Results representing predicted transverse tensile, longitudinal shear, and transverse shear stress-strain curves 
are presented, along with plots of the local fields indicating the damage progression within the microstructure. It is 
demonstrated that the damage model functions appropriately at the matrix scale, enabling localization of the damage 
to simulate failure of the composite material. The influence of the repeating unit cell geometry and the effect of the 
directionality of the applied loading are investigated and discussed. 


1. Introduction 

Progressive damage modeling of polymer matrix composites has received a great deal of attention in recent years 
as predictive capabilities for the complex nonlinear behavior of these materials are sought. A significant outstanding 
question is in regards to the proper, or most beneficial, scale on which the material damage model should function. 
Modeling the composite at the scale of the constituents (i.e., fiber and matrix materials) enables use of models for 
the matrix material that are initially isotropic and avoids ad-hoc coupling rules for the various damage components. 
However, a validated micromechanics method is required, and often, the constituent scale data best suited for charac- 
terizing the damage model are unavailable. In contrast, the macro-scale approach, wherein the composite is treated as 
an effective anisotropic material, will always be more computationally efficient, yet it requires a more complex damage 
model to handle the extreme anisotropy of the composite material with its directionally dependent multiple damage 
mechanisms. Further, the way in which the damage components interact (i.e., how much shear stiffness reduction 
results from a given amount of normal stiffness reduction) must be characterized, or, more commonly, assumed. 

Herein, an anisotropic, three-dimensional, progressive damage model, capable of modeling composites at the 
macro (composite) scale, as well as at the micro (constituent) scale, is presented. The development follows the 
approach outlined by Barbero [1] in the case of plane stress. The damage model uses the continuum damage approach 
first introduced by Kachanov [2, 3] in which continuous damage variables are introduced to quantify the stiffness 
reduction of a material in various directions. The theoretical background of this approach, and specifically, the elastic 
stiffness reduction damage approach, is discussed in detail by Lemaitre [5] and Chaboche [6], Kruch and Chaboche 
[18], Kajcinovic [7], Voyiadjis and Kattan [8], and Talreja [9, 10]. Further, Ladeveze [14, 15] and Pineda et al. [12] 
have reported progressive stiffness reduction damage models aimed at polymer matrix composites at the composite 
level. These models are motivated by the microstructure, but cannot consider micro scale effects like random fiber 
distribution and void content/morphology. Pineda et al. [11] and Bednarcyk et al. [13] presented stiffness reduction 
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damage models operating at the micro (fiber/marix) scale, but these studies did not examine micro scale architectural 
effects. Micromechanical investigations considering random fiber distributions of unidirectional (UD) fiber reinforced 
plastics (FRP) using an elasto-plastic isotropic damage model for the matrix material with the focus on fiber-matrix 
debonding were carried out by Melro et. al. [16, 17]. Chen and Ghosh [19] studied the strain rate dependent 
deformation and failure of composite materials at the micro-scale for different ideally packed unit cells under different 
strain rates and loading conditions. However, none of these works studied the influence of the coupling of the damage 
variables. 

The progressive damage model presented below is characterized for a polymer resin based on neat matrix ten- 
sile data. The model has been implemented within Abaqus [4] and applied to model the matrix material within a 
unidirectional polymer matrix composite. Three different composite repeating unit cells (with constant fiber volume 
fraction), representing three different fiber architectures, have been considered: hexagonal dense packing (i.e., equi- 
lateral hexagonal packing), and 5 and 15 fiber repeating unit cells (RUCs) with random fiber packing. The boundary 
conditions necessary to enforce periodicity at the RUC boundaries, while preserving the ability to apply arbitrary 
global (far-field) strains, are also presented. Results are shown for the three RUC microstructures for applied trans- 
verse tension, longitudinal shear, and transverse shear as a function of a damage model parameter that controls the 
extent of shear damage. The functionality of the damage model and its ability to predict damage localization, are 
highlighted, as are the effects of the directionality of the applied transverse tension and longitudinal shear loading, 
with respect to the RUC geometries. 


2. Progressive Damage Model 

A stiffness reduction progressive damage model has been developed and implemented based on the two-dimensional 
approach given by Barbero [1]. For application to the matrix material within a composite, a three-dimensional exten- 
sion of the model is needed due to the local three-dimensional stress state within a composite material. For a pure 
monolithic matrix, it is reasonable to consider the material to be initially isotropic, however when subjected to defor- 
mation, the formulation enables the material to exhibit anisotropic behavior depending on the activation (influence) of 
the damage mechanisms in the different directions. To allow for arbitrary damage, the in situ material model should 
have the capability to allow for different coupling of the damage variables, to represent, for example, penny shaped 
cracks or cavity type damage. The consequence is a very general formulation of the damage model, wherein the phys- 
ical interpretation of the material parameters is maintained. The model is based on the principle of energy equivalence 
[2] for small strains, and utilizes six scalar damage variables associated with the Young’s and shear moduli. Thus, the 
proposed Gibbs free energy is given by. 
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where p is the density of the material, oy, are the stress components, I) — [D\\D 22 D 33 D\ 2 D\ 3 D 23 \ are the dam- 
age variables, and E h G , s and v, ; are the pristine (undamaged) Young’s moduli, shear moduli, and Poisson’s ratios, 
respectively. The material compliance tensor, S, is then given by the second derivative of the Gibbs free energy with 
respect to the stress tensor. 
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The formulation relies on a damage surface, g(Y(D),y(<5)), and a damage potential, f(Y(D),y(6)), which are 
functions of thermodynamic forces, Y, and a hardening/softening function, y(&), with 6 being the damage harden- 
ing/softening variable. Moreover, it is associative, which requires / = g. When damage is accumulating. 


f =8=8-7 = 0 

with 

g = y/YHY 


(3) 

(4) 


7 = r + To = Ci(e c 2 - 1) + y 0 (5) 

where c\ and c 2 are material constants, y (l is the damage threshold, and H is a matrix of material parameters 
responsible for the interaction of the damage variables’ effect on the stiffness reduction of the material. If / is less 
than zero, damage does not evolve. 

The vector of thermodynamic driving forces is given by 


and the damage evolution is 


v_ 

P dD 




(6) 


(7) 


where the Lagrangian (damage) multiplier A can be interpreted as the damage rate magnitude and is the direc- 
tion of damage evolution in y-space. Further, the damage rate magnitude affects the hardening variable rate as 


Making use of the consistency condition 
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the Lagrange (damage) multiplier can be obtained 
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where e is the strain tensor. Then, the damage variable evolution can be calculated from (7). Finally, from (3), (5) 
and (8), the hardening variable evolution is 


6 = -A (11) 

A standard return mapping algorithm with an elastic predictor/damage corrector strategy was used for the numer- 
ical implementation into the UMAT user material routine of the commercial FE software Abaqus [4]. 

3. Damage Interaction 

Based on (7) and (3-5), the evolution of one damage component is 

bi = = ^[2 H n Y, + (H ]2 H 2] )Y 2 + (H n H v )Y 3 + {H u H Al )Y A + (H l5 H 5l )Y 5 + (H l6 H 6l )Y 6 \ (12) 

dY 2 g 

Thus, a uniaxial stress, which will produce only one component of Y (see (6) and (1)), denoted by Yj, will induce 
damage in 1 -direction, provided Hij is non-zero. Thus, to simulate a fracture mechanics type interaction (penny 
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shaped crack), Hu, Hu, and Hu can be chosen as non-zero. Then, crn, cr^, and cr 13, would induce damage that 
would cause a stiffness reduction in the 11-, 12- and 13-components. Conversely, if the damage evolution of the 
material is assumed to be isotropic, then the components H t j can be chosen as. 


H n = H 2 2 = W 33 = H nn 

Hu = ftu = H 2A = H 42 = H 15 = H S1 = H 35 = H 53 = H 26 = H 62 = H 36 = H 63 = H sn (13) 

Hu — H 33 — i/66 — H ss 

In addition to specification of these H matrix terms, c 1, c 2 and y (l must be chosen to fully characterize the damage 
model. 

4. Damage Model Characterization 

The damage model was characterized for tension through correlation with a cyclic tension test on a neat epoxy 
resin material (RIMR 135/RIMH 137, mass mixing ratio 10:3, 8 hr curing time at 40° C). The epoxy was vacuum 
injected into a standard dogbone mold to minimize voids. Cyclic tension was achieved using a constant crosshead 
speed of 10 mm/s, with strains measured via a digital image correlation system. 



0 0.5 1 1.5 2 2.5 3 3.5 4 

Strain [%] 

Figure 1 : Model correlation with epoxy tensile test data 


The model characterization is shown in Figure 1. Note that the test specimen failed suddenly during the final 
reloading cycle at the point indicated. Furthermore, the nonlinear resin response appears to be due to both inelastic 
deformation and damage mechanisms as some permanent deformation is observed after unloading, albeit small. In the 
present formulation all nonlinear behavior is assumed to be associated with damage only. Since the damage model is 
based purely on stiffness reduction, upon unloading, there is zero permanent strain. This is a reasonable approximation 
of the actual epoxy tensile response as only a small amount of permanent strain was evident upon unloading. Based 
on the correlation shown, the isotropic elastic properties and damage model parameters obtained for the epoxy are 
shown in Table 1 . 

The Poisson’s ratio of the epoxy given in Table 1 was not measured, but set to a typical value. The influence of the 
shear damage parameter, H ss , will be examined via parametric studies on a unidirectional composite repeating unit 
cell. 
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Table 1 : Material parameters for RIMR 135/RIMH 137 epoxy matrix 


Matrix 

E 

3000 MPa 

V 

0.38 

Cl 

27.5 MPa 

Cl 

-2.8 

Hfin 

0.75 

To 

0.45 MPa 


5. Finite Element Models 

The three-dimensional damage model presented above has been implemented into Abaqus [4] as a UMAT user- 
defined material constitutive routine. To represent a unidirectional composite material, three RUC geometries have 
been discretized using Abaqus, as shown in Figure 2. 





Figure 2: Abaqus finite element model RUCs used to represent a 0.49 fiber volume fraction unidirectional composite. Left: HDP, center: 5 fiber 
random, right: 15 fiber random 


Figure 2 left shows a microstructure with hexagonal dense packing (HDP) meshed with 2,166 C3D8 elements. 
Figure 2 center contains 5 fibers with random locations meshed with 4,420 C3D8 elements, and Figure 2 right contains 
15 fibers with random locations meshed with 13,202 C3D8 elements. In all cases, the fiber volume fraction is 0.49. 
The fiber is assumed to be carbon; the transversely isotropic elastic properties are shown in Table 2 


Table 2: Material parameters for carbon fiber; 1 -direction is longitudinal fiber direction 


Fiber 

Ei 

290 GPa 

Ei 

20 GPa 

Vl2 

0.2 

Gn 

20 GPa 

G 13 

9 GPa 


The influence of external boundaries is removed by imposing periodic boundary conditions. Within Abaqus, this 
is accomplished using constraint equations, as follows for a parallelepiped model. First, node sets for every node 
on each face, edge, and vertex must be generated. The nodes at the vertices must not be included in the node sets 
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associated with edges or faces. Similarly, the nodes at the edges are excluded from the face node sets. Next, the node 
pairs (two nodes that are coupled via a constraint equation) for the three sets of faces, and the node groups for the 
three sets of edges, must be identified. The node pair for a face node is always located on the opposite face. For 
example, referring to Figure 3, a node on face 4 with the coordinates (x = xi,y = 0,z = Zi) has its partner node on 
face 2 with (x = x\,y = b, z = Zi ), where b is the length of the RUC in the v-dircction. Similarly, the node groups 
for edges must be determined. The edges are categorized based on direction; for instance, the x -direction edge set 
includes edges 9, 10, 11 and 12. Then the node groups within each set of edges are found at the same coordinate along 
the edge direction. For example, a node on edge 9, with the coordinates (x — xi , y = 0, z = 0) is grouped with nodes 
on edges 10 (x = x\ ,y = b,z = 0), 11 (x = x\ ,y — b,z — c) and 12 (x = x\,y - 0,z = c), where c is the length of the 
RUC in the z-dircction. 


b 



Z 



Z 



Figure 3: Reference numbering for application of periodic boundary conditions on a parallelepiped finite element model using equation constraints 
within Abaqus. Left: faces, center: vertices, right: edges 


In Abaqus, when coupling nodes via equation constraints, it is only possible to couple the displacements (as 
opposed to forces). Using three reference points (RPs), designated as RPX at the position of vertex 1, RPY at the 
position of vertex 6, and RPZ at the position of vertex 8, the complete (arbitrary) far field strain state can be specified 
by applying displacement boundary conditions to these three RPs. 

The equation constraints must be defined such that, for every node pair/group, only one degree of freedom remains 
in the system, without ever eliminating the same degree of freedom more than once. Denoting the displacement in 
coordinate direction i and vertex, edge, face, or reference point j as u J , one possibility that accomplishes this is 
For edges: 
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For faces: 
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For vertices: 
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where these equation constraints are imposed on paired and grouped nodes, as described above. Thus, in order 
to apply a far-held strain, only the displacements of the reference points (RPs), according to the strain state, must be 
specified as boundary conditions. Figure 4 shows several examples of deformed RUC plots, where the above periodic 
boundary condition application strategy has been employed. As can be seen, the approach maintains periodicity, even 
in the presence of significant damage localization. 




Figure 4: Example deformed plots (deformation scale factor = 5) showing that periodicity is maintained, even in the presence of significant damage 
localization 



It is well-known that, when applying damage models such as this in finite element analysis, pathological mesh- 
dependence will be present once localization occurs [21, 22]. Thus, for each RUC geometry, it was desired to keep 
the element size constant during the study. Further, the global element size variation was kept small by restricting 
the element volume during the discretization process in Abaqus. The element volume and was set to the same value 
as used in the material calibration step. In the future, length scaling will be added to the model to ensure consistent 
energy dissipation in the presence of varying mesh sizes [24]. 
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6. Results and Discussion 


The impact of the RUC geometry (see Figure 2) on the predicted local and global response of the composite 
has been examined while also varying the value of the shear damage parameter H ss . For the cases examined, the 
shear-normal interaction damage parameters, H m , has been set to zero, eliminating local shear stiffness reduction due 
to normal stresses as well as normal stress reduction due to shear stresses. Figure 5 shows the predicted transverse 
tensile response of the composite as a function of the RUC geometry (see Figure 2) and the direction (x or y) of the 
applied loading for the case of H ss — 5. Maximum principal strain fields, which highlight the location of damage 
within the RUC, are plotted in Figure 6 for five values of the applied global strain. Note that these global strain values 
are not exact since adaptive time steps were employed in the finite element simulations. 



Figure 5: Comparison of the predicted effective composite transverse tensile response as a function of RUC geometry and direction of loading 
(solid line = loading applied in x-direction, dashed line = loading in the y direction) for the case of H ss = 5 

While all cases predict nearly identical initial elastic stiffness, once damage initiates, there are significant dif- 
ferences in the global deformation (i.e., damage) response due to microstructural variation in the three RUC repre- 
sentations. For the HDP RUC, there is a large ( 22%) difference between tensile loading in the x-direction and the 
y-direction in terms of the stress at which departure from linearity occurs and at which a stress drop occurs. This stress 
drop is associated with the onset of damage localization within the RUC. As evident in Figure 6, the local fields are 
different (even initially), which leads to initiation of damage localization much earlier when the HDP RUC is loaded 
in the x-direction. This is in spite of the fact that, elastically, the RUC behaves transversely isotropically, implying 
that the average elastic fields are identical for tension in either direction. The final damage patterns are also quite 
different, with localization progressing completely through the RUC much earlier in the case of y-direction loading as 
compared to loading in the x-direction. 

For the 5 fiber random RUC, Figure 5 shows that, while the difference in the x- and y-direction tensile curve is 
reduced compared to the HDP RUC, it is still significant ( 18%). The y-direction loading results in earlier onset of 
damage localization because of a concentration that arises in the matrix between two fibers that are close together in 
the y-direction, as can be seen in Figure 6. Both 5 fiber RUC loading simulations exhibit earlier damage localization 
compared to the HDP RUC simulations due to such concentrations between nearby fibers. In contrast, both 15 fiber 
RUC loading simulations predict localization at nearly the same time (strain). The damage patterns shown in Figure 
6 are qualitatively similar for each level of applied global strain, with the x-direction loading exhibiting a somewhat 
greater extent of damage localization at the applied global strain of 1%. This is manifested in the somewhat more 
compliant response after damage initiation in Figure 5. The fact that the x- and y-direction responses are becoming 
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closer indicates that, as one would expect, as a larger RUC with more fibers is considered, the model provides a better 
representation of the actual (transversely isotropic) material response. 

Figure 7 shows the predicted x-direction tensile response of the 15 fiber RUC as a function of the value of the 
shear damage parameter, H ss , while the corresponding maximum principal strain fields are plotted (at an applied 
global strain level of approximately 1%) in Figure 8. It is clear from both figures that, by increasing H ss , which 
causes an increase in the shear damage induced by local shear stresses, the progression of damage is accelerated in 
the composite. This highlights the importance of capturing normal-shear coupling in the micromechanics solution in 
the presence of damage and the value of shear test data for the neat resin. 

Figure 9 and Figure 10 show the response and local maximum principal strain plots for applied global longitudinal 
shear strains, y xz and y y7 . All cases show good agreement in the initial linear regime, with the 5 fiber RUC loading in 
yz shear exhibiting slightly higher stiffness. 

In general, the effect of the orientation of the applied shear strain (xz vs. yz) has much less of an effect on 
the longitudinal shear response compared to the effect of the orientation of the tensile loading (x vs. y) on the tensile 
response. The effect is, however, still present, and, once again, it is minimal in the case of the 15 fiber RUC, suggesting 
the need to model a good number of fibers within an RUC (to approach a converged representative volume element 
(RVE)). For the HDP RUC, while localization initiates at roughly the same time for both xz and yz shear, the damage 
progresses much more quickly under yz shear loading. It is clear from the damage localization patterns in Figure 10 
that this is due to the more tortuous path of the damage in the case of xz loading, which leads to the tougher response 
in Figure 9. In the 5 fiber RUC, the opposite trend is present, with the damage progressing more rapidly, with a 
straighter path, under xz shear loading. The effect of the RUC geometry itself on the longitudinal shear response is 
also still present, with nearby fibers leading to earlier damage initiation for the 5 fiber and 15 fiber RUCs. However, 
the RUC effect is muted compared to the transverse tensile response. 

Figure 11 and Figure 12 show the effect of the matrix shear damage parameter, H ss , on the predicted xz shear 
response of the 15 fiber RUC. For H ss = 0, no shear damage can occur, and the effective response remains linear to 
very high strains when damage due to local normal strains can arise. 

As one would expect, increasing H ss , which increases the rate of local matrix shear stiffness reduction in response 
to local shear stress, has a primary effect on the predicted composite longitudinal shear response. The effect is much 
greater than the effect of H ss on the predicted transverse tensile response (see Figure 7). Figure 13 and Figure 14 show 
the predicted transverse (xy) shear response of the composite, along with plots of the local maximum principal strain 
fields. The 5 fiber RUC predicts slightly higher elastic transverse shear stiffness compared to the HDP and 15 fiber 
RUCs, which are in close agreement in the linear region. The magnitude of the effect of the RUC geometry on the 
predicted transverse shear stress-strain response of the composite is similar to its effect on the predicted longitudinal 
shear stress-strain response. As before, local stress concentrations arising in the matrix between nearby fibers tend 
to cause early onset of localization due to damage. Finally, all three RVE geometries give rise to straight damage 
localization paths, leading to similar post-peak xy stress-strain responses. 

7. Conclusion 

A three-dimensional, anisotropic, progressive damage model, relying on six scalar damage variables, has been 
developed following the approach outlined by Barbero [1] in the case of plane stress problems. By extending the 
model to three-dimensional stress states, it is now capable of modeling the behavior of the constituents within a 
composite, wherein the stress state is inherently three-dimensional. Further, because it is anisotropic, the model could 
also be applied to any effective anisotropic material; for example the tows within woven composites (though not done 
herein), which also experience three-dimensional stress states. The damage model is also arbitrary in terms of the 
interaction among the six damage variables and the stress components, which is controlled by the chosen form of 
the damage potential. The damage model has been implemented within Abaqus [4] as a user-defined constitutive 
model. The model was characterized for a neat resin material and then applied to model the transverse tensile, 
longitudinal shear, and transverse shear response of a unidirectional composite composed of that matrix and carbon 
fibers. Finite element models representing RUC geometries with HDP, 5 randomly placed fibers, and 15 randomly 
placed fibers were considered. Periodicity conditions were imposed on the boundaries of the finite element RUC 
model to ensure representation of the composite material. Results showed that the damage model is capable of 
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predicting damage localization within a periodic composite in response to both normal and shear loadings for all 
three RUCs. Further, while the linear portion of the predicted stress-strain curves of all three repeating cells were 
nearly coincident, significant differences were present once damage initiated due to variations in the represented 
microstructures. The HDP and 5 fiber RUCs also exhibited significant differences based on the directionality of the 
applied transverse tension and longitudinal shear loading. In contrast, the 15 fiber RUC showed only minor influence 
of the loading direction. The influence of the shear damage parameter on the predictions of the 15 fiber RUC was also 
investigated and shown to significantly soften the post damage initiation response for both applied transverse normal 
and longitudinal shear loading. 
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Figure 6: Predicted local maximum principal strain plots during simulated transverse tensile tests at different applied global strains e as a function 
of RUC geometry and direction of loading for the case of H ss = 5. The deformation scale factor is 5 
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Figure 7: Comparison of the predicted effective composite X-direction transverse tensile response for the 15 fiber RUC as a function of the shear 
damage parameter, H ss 
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Figure 8: Predicted local maximum principal strain plots during simulated x-direction transverse tensile tests on the 15 fiber RUC as a function of 
the shear damage parameter, H ss , at an approximate applied global strain of e = 1%. The deformation scale factor is 5 
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Figure 9: Comparison of the predicted effective composite longitudinal shear response as a function of RUC geometry and direction of loading for 
the case of H ss = 5 
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Figure 10: Predicted local maximum principal strain plots during simulated longitudinal shear tests at different applied global shear strains y as a 
function of RUC geometry and direction of loading for the case of H ss = 5. The deformation scale factor is 5 
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Figure 11: Comparison of the predicted effective composite xz longitudinal response for the 15 fiber RUC as a function of the shear damage 
parameter, H ss 



Figure 12: Predicted local maximum principal strain plots during simulated xz longitudinal shear tests on the 15 fiber RUC as a function of the 
shear damage parameter, H ss , at an approximate applied global shear strain of gamma = 1.35%. The deformation scale factor is 5 
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Figure 13: Comparison of the predicted effective composite transverse (xy) shear response as a function of RUC geometry for the case of H ss = 5 
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Figure 14: Predicted local maximum principal strain plots during simulated transverse (xy) shear tests at different applied global shear strains y as 
a function of RUC geometry and direction of loading for the case of H ss = 5. The deformation scale factor is 5 
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